function [ outdata ] = ratheartsegmentation( basedir, subdirs )
% ratheartsegmentation: segment myocardium in 6 part; draw graphs over time
%
% Myocardium segmentation of the left ventricle for short axis MR scans.
% Generates a curve showing the evolution of R2 values within each region
% over time.
%
% Usage: 
%    outdata = ratheartsegmentation(basedir, subdirs);
%
% where:
%    basedir = the dicom directory of a specific rat
%    subdirs = a cell array of strings with subdirectories pointing to the
%              T2 subdirectories.
%    outdata = a struct array, containing for each day some relevant data.
%              outdata(day).Rmap contains the registered Rmap for
%              example.

% Adapted by Stefan Klein. Based on code by Gabrielle Giorico, August 2011.


%Original image

%path = 'C:\Users\ggiorico\Documents\Segmentation\OriginalImages\';
path = 'D:\data\proj\ratheart\repro\20110815\1092\T2\012 FSE-TE 26.4\';
testdaypath = 'D:\ck\matlab\gabrielle\outdir2\';
%info = dicominfo(strcat(path,'Rat', int2str(nb_rat),'\fixed.dcm'));

% determine postproc dirs.
nb_days = numel(subdirs);
for day = 1:nb_days
  studydir = fullfile(basedir, subdirs{day} );
  % remove slashes at the end
  while ( studydir(end) == '/' || studydir(end) == '\' )
    studydir = studydir(1:(end-1));
  end
  outdir{day} = [ studydir 'postproc\' ];
  if ~exist(outdir{day}, 'dir')
    error( 'ratheartsegmentation:invalidoutputdir', '%s : invalid postproc dir. Run first ratheartT2fit!', outdir{day});  
  end
end

% all results are written to the first day's postproc directory
resultdir = outdir{1};

% load the first two clicked points from the first day
points = load( fullfile( resultdir, 'points.mat') );
points = points.points( [1 2], : );
% store x and y coordinates separately, to use gabrielle's code.
get_X = points(:,1)';
get_Y = points(:,2)';

% Load the first day's R map
cd(resultdir)
d = dir;
for j=3:size(d,1)
    x=d(j).name;
    if size(char(regexp(x,'R1map', 'match')),2)>=1 %read in the fixed longest IR image
        info_zeroline = dicominfo( fullfile( resultdir, 'R1map.dcm' ) );
    elseif size(char(regexp(x,'R2map', 'match')),2)>=1 %read in the fixed longest IR image
        info_zeroline = dicominfo( fullfile( resultdir, 'R2map.dcm' ) );
    end   
end
original_map_zeroline = double(dicomread(info_zeroline));

% compute the mask from the R map and the centre of mass of the mask
mask = double( original_map_zeroline ~= 0 );
S = sum( mask(:) );
[r,c]= ndgrid(1:size(original_map_zeroline,1),1:size(original_map_zeroline,2)); %array of coordinates
Xo= ( sum(mask(:).* c(:)) / S );
Yo= ( sum(mask(:).* r(:)) / S );

% coordinates of the middle of the segment formed by points 1 and 2.
m_X = (get_X(1,1)+ get_X(1,2))/2;
m_Y = (get_Y(1,1)+ get_Y(1,2))/2;

% compute angle of zero line
m_angle = atan2(( m_Y - Yo), (m_X - Xo));
if(m_angle<0)
    m_angle = m_angle + 2*pi;
end
m_angled = (180 * (m_angle))/pi;

% Initialization
array_meanr = []; 
array_std = [];
array_crlb = [];

% Loop over all the days
for day=1:nb_days
        
    cd(outdir{day})
    d = dir;
    for j=3:size(d,1)
        x=d(j).name;
        if day ==1
            if size(char(regexp(x,'R1map', 'match')),2)>=1 %read in the fixed longest IR image
                infoRmap = dicominfo( fullfile( outdir{day}, 'R1map.dcm') );
                infoRerrormap = dicominfo( fullfile( outdir{day}, 'R1errormap.dcm') );
                infoTE = dicominfo( fullfile( outdir{day}, 'FirstIm.dcm') );
            elseif size(char(regexp(x,'R2map', 'match')),2)>=1 %read in the fixed longest IR image
                infoRmap = dicominfo( fullfile( outdir{day}, 'R2map.dcm') );
                infoRerrormap = dicominfo( fullfile( outdir{day}, 'R2errormap.dcm') );
                infoTE = dicominfo( fullfile( outdir{day}, 'FirstIm.dcm') );
            end
        else
             if size(char(regexp(x,'registeredR1map', 'match')),2)>=1 %read in the fixed longest IR image
                infoRmap = dicominfo( fullfile( outdir{day}, 'registeredR1map.dcm') );
                infoRerrormap = dicominfo( fullfile( outdir{day}, 'registeredR1errormap.dcm') );
                infoTE = dicominfo( fullfile( outdir{day}, 'registeredFirstIm.dcm') );
            elseif size(char(regexp(x,'registeredR2map', 'match')),2)>=1 %read in the fixed longest IR image
                infoRmap = dicominfo( fullfile( outdir{day}, 'registeredR2map.dcm') );
                infoRerrormap = dicominfo( fullfile( outdir{day}, 'registeredR2errormap.dcm') );
                infoTE = dicominfo( fullfile( outdir{day}, 'registeredFirstIm.dcm') );
             end
        end   
    end
    
    Rmap= double(dicomread(infoRmap));
    Rerrormap= double(dicomread(infoRerrormap));
    TE = double(dicomread(infoTE));
    mask = int16(Rmap~=0);

    % center of gravity of the binary mask; SK: use one from the first day.
%     MaskS=sum(mask(:));
    %[r,c]= ndgrid(1:size(mask,1),1:size(mask,2)); %array of coordinates
%     Xc= ( sum(Mask(:).* c(:)) / MaskS );
%     Yc= ( sum(Mask(:).* r(:)) / MaskS );

    Xc = Xo;
    Yc = Yo;

    %coordinates of non-zeros pixels 
    %ArrayofAngle=[];
    
    for y=1:size(mask,1)
        for x=1:size(mask,2)
            if(mask(y,x)~=0)

                angle = atan2( (y-Yc), (x-Xc) ); %result in radian
                if(angle<0)
                    angle = angle + 2*pi;
                end
                angled = (180 * (angle))/pi;  %result in degrees [0,360]

                % Regions for the segmentation of the basal and mid short axis
                % For the apical short axis, only 4 regions
                mask(y,x)= floor(mod(angled + 360 - m_angled, 360) / 60) + 1;
            end
        end
    end

    % To find the values corresponding to a region in the Rmap
    vec_region_day = [];
    vec_std_day = [];
    vec_crlb_day = [];

    for r=1:6  %region 1 to 6
        rvalues = Rmap( (mask == r) & (Rmap>0) ) ;
        %rvalues = rvalues(rvalues>0); % get rid of negative values.
        rvalues_final = rvalues/10000;
        rerrvalues = Rerrormap( (mask == r) & (Rmap>0) );
        rerrvalues_final = rerrvalues/10000;

        % store mean r, std r, and mean crlb
        vec_region_day = [vec_region_day, mean(rvalues_final)];
        vec_std_day = [vec_std_day, std(1./rvalues_final)];
        vec_crlb_day = [vec_crlb_day, mean(rerrvalues_final)]; 
       
    end % end of the region loop
    
    % Update the matrix containing the mean value of the region per day
    % (each line is a day)
    array_meanr = [array_meanr; vec_region_day];
    array_std = [array_std; vec_std_day];
    array_crlb = [array_crlb; vec_crlb_day];
    
    % store for user:
    outdata(day).Rmap = Rmap;
    outdata(day).Rerrormap = Rerrormap;
    outdata(day).firstTE = TE;
    outdata(day).postprocdir = outdir{day};
    outdata(day).segmentation = mask;
    
end % End of the day loop

% Plot
screensize = get(0,'Screensize');
figure('Name','Analysis results','Position', [0.02*screensize(3),0.1*screensize(4),0.95*screensize(3),0.5*screensize(4)]);
vec_xaxis = [];

% TO DO: Look for the command lign that create the vector directly
for i=1:nb_days
    vec_xaxis = [ vec_xaxis, i];
end

% Graphs
shiftbetweenstdandcrlb = 0.05;
% ylim = [min(min(1./array_meanr))-2*max(max(array_std)) max(max(array_meanr))+2*max(max(array_std))];
ylim = [0 100];% 1.9*max(max(1./array_meanr))];
regionnames = {'anteroseptal','anterior','anterolateral','inferolateral','inferior','inferoseptal'};
for r=1:6
    subplot(2,3,r);
    errorbar( vec_xaxis, 1./array_meanr(:,r)', array_std(:,r), '.b-'); 
    hold on;
    if r==1
        % dummy, to fake legend
        errorbar( 10, 1, 1, '.r-'); 
    end
    errorbar( vec_xaxis + shiftbetweenstdandcrlb, 1./array_meanr(:,r)', (array_crlb(:,r)./(array_meanr(:,r).^2))', '.r'); 
    axis([ 0.8 nb_days+0.2 ylim]);
%     h = title(['Region' num2str(r)]);
      h = title(regionnames(r));
    pos = get(h,'Position');
    pos(2) = ylim(1)+ (ylim(2)-ylim(1))*0.85;
    set( h, 'Position', pos);
    xlabel('Days');
    set(gca,'XTick', 1:nb_days );
    ylabel('T2 values [ms]');
    if r==2        
        lh=legend({'st.dev.', 'CRLB'}, 'Location', 'North','Orientation', 'horizontal');
    end
end
pos = get(lh, 'Position');
pos([1 2]) = [0.3 0];
set(lh, 'Position', pos);

% 1./array_meanr(:,2)'



figure('Name', 'Region definition', 'position', [0.05*screensize(3),0.6*screensize(4),0.25*screensize(3),0.3*screensize(4)]); 
imshow( outdata(1).segmentation, [],'InitialMagnification','fit'); 
colorbar;
title( 'Region definition at first day.' );

figure('Name', 'Inspect registered images', 'position', [0.3*screensize(3),0.45*screensize(4),0.63*screensize(3),0.42*screensize(4)]);
for day = 1:nb_days
    subplot(2,nb_days, day);
    imshow( outdata(day).firstTE, [],'InitialMagnification','fit' );
    title( ['first TE, day ' num2str(day) ] );
    subplot(2,nb_days, day + nb_days);
    imshow( outdata(day).Rmap, [],'InitialMagnification','fit' )
    title( ['Rmap, day ' num2str(day) ] );
end

outputfile=[resultdir 'Results_Repro.txt'];
fid = fopen(outputfile, 'a');
% fprintf(fid, 'T2_Spio T2CR_Spio T2STD_Spio T2_Myocard T2CR_Myocard T2STD_Myocard R2_Spio R2CR_Spio R2STD_Spio R2_Myocard R2CR_Myocard R2STD_Myocard\n');
fprintf(fid, 'Day T2_r1 T2_r2 T2_r3 T2_r4 T2_r5 T2_r6 T2err_r1 T2err_r2 T2err_r3 T2err_r4 T2err_r5 T2err_r6 R2_r1 R2_r2 R2_r3 R2_r4 R2_r5 R2_r6 R2err_r1 R2err_r2 R2err_r3 R2err_r4 R2err_r5 R2err_r6 \n');
for day=1:nb_days
    fprintf(fid, '%i %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n', ...
    day,...
    1/(array_meanr(day,1)),...
    1/(array_meanr(day,2)),...
    1/(array_meanr(day,3)),...
    1/(array_meanr(day,4)),...
    1/(array_meanr(day,5)),...
    1/(array_meanr(day,6)),...
    array_crlb(day,1)/(array_meanr(day,1)^2),...
    array_crlb(day,2)/(array_meanr(day,2)^2),...
    array_crlb(day,3)/(array_meanr(day,3)^2),...
    array_crlb(day,4)/(array_meanr(day,4)^2),...
    array_crlb(day,5)/(array_meanr(day,5)^2),...
    array_crlb(day,6)/(array_meanr(day,6)^2),...
    (array_meanr(day,1)),...
    (array_meanr(day,2)),...
    (array_meanr(day,3)),...
    (array_meanr(day,4)),...
    (array_meanr(day,5)),...
    (array_meanr(day,6)),...
    array_crlb(day,1),...
    array_crlb(day,2),...
    array_crlb(day,3),...
    array_crlb(day,4),...
    array_crlb(day,5),...
    array_crlb(day,6));
end
    disp( ['Saved results to ' outdir{1} 'Results_Repro.txt'] );
    fclose(fid);


% subplot(2,3,2);
% errorbar(vec_xaxis, array_meanr2(:,2)', array_std(:,2),'g');legend('Region2');
% axis([ 0.8 size(array_meanr2(:,2)',2)+0.2 0 0.15 ]);
% xlabel('Days');
% ylabel('R2 values [ms]');
% subplot(2,3,3);
% errorbar(vec_xaxis, array_meanr2(:,3)', array_std(:,3),'r');legend('Region3');
% axis([ 0.8 size(array_meanr2(:,3)',2)+0.2 0 0.15 ]);
% xlabel('Days');
% ylabel('R2 values [ms]');
% subplot(2,3,4);
% errorbar(vec_xaxis, array_meanr2(:,4)', array_std(:,4),'m');legend('Region4');
% axis([ 0.8 size(array_meanr2(:,4)',2)+0.2 0  0.15 ]);
% xlabel('Days');
% ylabel('R2 values [ms]');
% subplot(2,3,5);
% errorbar(vec_xaxis, array_meanr2(:,5)', array_std(:,5),'k'); legend('Region5');
% axis([ 0.8 size(array_meanr2(:,5)',2)+0.2 0 0.15 ]);
% xlabel('Days');
% ylabel('R2 values [ms]');
% subplot(2,3,6);
% errorbar(vec_xaxis, array_meanr2(:,6)', array_std(:,6),'c'); legend('Region6');
% axis([ 0.8 size(array_meanr2(:,6)',2)+0.2 0 0.15 ]);
% xlabel('Days');
% ylabel('R2 values [ms]');


   